Coulomb crystallization in expanding laser-cooled neutral plasmas 
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We present long-time simulations of expanding ultracold neutral plasmas, including a full treat- 
ment of the strongly coupled ion dynamics. Thereby, the relaxation dynamics of the expanding 
laser-cooled plasma is studied, taking into account elastic as well as inelastic collisions. It is demon- 
strated that, depending on the initial conditions, the ionic component of the plasma may exhibit 
short-range order or even a superimposed long-range order resulting in concentric ion shells. In 
contrast to ionic plasmas confined in traps, the shell structures are built up from the center of the 
plasma cloud rather than from the periphery. 

PACS numbers: 52.27.Gr,32.80.Pj,52.38.-r,34.60.+z 



It is well known that, depending on the Coulomb cou- 
pling parameter (CCP) T = e 2 /ak B T, a plasma may 
show long-range order, short-range order or no order at 
all. In general, ordering effects can be expected in the 
so-called strongly coupled regime (r> 1) where the in- 
terparticle Coulomb interaction e 2 /a dominates the ther- 
mal energy k B T of the plasma particles. This parameter 
regime has been studied extensively in nonneutral plas- 
mas of laser-cooled ions confined in ion traps 
On the other hand, much less is known about the dynam- 
ics of finite, strongly coupled neutral plasmas without 
confinement. Due to their expansion, these plasmas are 
in a non-equilibrium state at all times, and it is not clear 
whether dramatic ordering effects such as Coulomb crys- 
tallization known from static trapped ionic plasmas can 
be observed in such a system. 

Experimentally, cold neutral plasmas could be real- 
ized only recently by photoionizing a cloud of ultracold 
(T <C IK) atoms Yet, under the present exper- 

imental conditions, the regime of a strongly correlated 
plasma cannot be reached jfj, 111! . Since the plasma 
is created in a completely uncorrelated state, the subse- 
quent conversion of potential into kinetic energy rapidly 
heats both the electron and ion subsystem, suppressing 
the development of substantial correlations |sl llCj. Ad- 
ditionally, inelastic collisions with Rydberg atoms, previ- 
ously formed b y th ree-body recombination (TBR), and 
the TBR itself |l2| heat the electron gas. Therefore, 
Tj decreases to unity while T e becomes even smaller, 
preventing strong correlation effects. However, further 
Dopplcr cooling of the ions during the plasma expansion 
has been suggested as a possible route to strongly coupled 
ultracold plasmas [lj, |l4j . 

In the following we provide the first theoretical de- 
scription for the expansion of such a laser-cooled neutral 
plasma. As we will demonstrate, strong-coupling phe- 
nomena can indeed occur in such a system, leading to 
the formation of surprisingly differentiated patterns in 
the ionic plasma given the appropriate initial conditions. 
From a nonlinear dynamics point of view one might char- 
acterize these phenomena as self organization of a system 



in a non-equilibrium state. 

In a first step we describe the collisionless plasma dy- 
namics by a set of coupled Vlasov equations for the elec- 
trons and ions neglecting any correlation effects. 
Cooling of the ions is modelled by introdu cing a Fokker- 
Planck term into the ion kinetic equation [16( 
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is the one-particle distribution function of the 
is the ion mass and the damping rate and 
the Doppler temperature T c are determined by the de- 
tails of the laser cooling process. Assuming quasineu- 
trality together with an adiabatic treatment of the elec- 
trons, one can show that the resulting ionic kinetic equa- 
tion permits a Gaussian selfsimilar solution /, (r, v) cx 
cxp(— r 2 /2a 2 — mi (v — 71*) /2k riTi ) l also found in the 
free plasma expansion problem [ill ITil Ht| . The rms- 
radius a of the plasma cloud, the hydrodynamic velocity 
parameter 7 and the temperatures Tj and T e evolve ac- 
cording to 
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As seen from Eq. (J2J), the action of the cooling laser is 
twofold. The ion temperature is driven towards its equi- 
librium value T c , while 7 is linearly damped out on a 
timescale of Eq. (0) still has one integral of motion, 
a 2 T e = const., reflecting the adiabatic electron cooling 
during the plasma expansion. If T» <C T e , an asymptotic 
description of the plasma dynamics can be obtained by 
neglecting terms of order 7 compared to the damping 
rate f3 in Eq. (J2J, which yields in the long-time limit 
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This is in marked contrast to the free expansion, where 
a 2 w cr 2 (0) + kB ll [a) t 2 'k°° t 2 . Therefore, continued 
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FIG. 1: The size of a laser-cooled plasma as compared to that 
of a freely expanding plasma of the same initial-state param- 
eters (dotted). The solution of Eq. (|5J (dashed) matches the 
PIC-treatment (circles), showing that the difference to the 
hybrid-MD simulation (solid) comes from the ionic correla- 
tions. The inset shows a comparison of the numerical solution 
with the analytical approximation Eq. (J^J. 



Dopplcr cooling of the ions not only reduces their temper- 
ature but also drastically retards the decrease in plasma 
density, supporting development of strong correlations. 

In a second step, we have performed a more elaborate 
numerical simulation based on a hybrid method treating 
the two plasma components on different levels of sophis- 
tication. Since the electrons are not strongly coupled and 
their relaxation time, determined by the electron plasma 
frequency, is small compared to both the plasma expan- 
sion time and the inverse of the ionic plasma frequency 
Wp = i k /47re 2 p/mi, an equilibrium fluid model provides 
an adequate description of the electron dynamics 
We account for the initial electron evaporation by deter- 
mining the fraction of trapped electrons from the results 
of Ref. 0. The ions move as classical particles under the 
action of the electronic mean field and the direct ion-ion 
interaction calculated with a particle tree-procedure de- 
veloped in |l8|. In analogy to Eq. cooling is described 
by adding a Langevin force to the ion equations of mo- 
tion. The electron temperature is determined by energy 
conservation for the total system consisting of the plasma 
and the radiation field. This hybrid treatment allows us 
to study effects of strong ion correlation over long times, 
since atomic timescales need not be resolved as in a full 
molecular dynamics (MD) simulation |9Ul9|. 

In Fig. n we compare the time dependence of the 
plasma rms-radius obtained from Eq. J5J with the hybrid- 
MD simulation and the analytical approximation Eq. J3J). 
There is good overall agreement between the two numer- 
ical approaches. Moreover, they both nicely reproduce 
the asymptotic %/i-behavior of Eq. ©. On the other 
hand, the width calculated from the MD simulation is 
significantly shifted to lower values. A comparison with 
a particle-in-cell (PIC) treatment of the ions, also shown 



in Fig.Q] clearly reveals that the slower plasma expansion 
is due to the negative correlation pressure 0, 0] , which 
partly compensates the thermal electron pressure. Note 
that here the influence of ion correlations is completely 
different from the case of free plasma expansion, where 
the initial correlation heating was found to dominate the 
negative correlation pressure and hence accelerates the 
plasma expansion |17j . 

Up to this point, we have taken into account the 
electron-ion interaction on the basis of a mean field de- 
scription only. However, it has been found that electron- 
ion collisions leading to the formation of Rydberg atoms 
through TBR may significantly alter the expansion dy- 
namics at these low electron temperatures [y, [llj . In or- 
der to include these processes in our description, we use a 
Monte-Carlo treatment 0, H(| to account for TBR and 
inelastic electron-Rydberg atom collisions. In addition 
to these processes the influence of the cooling laser on 
the Rydberg atoms should be addressed. By the very 
nature of the cooling process, a significant fraction of 
the ions is found in an excited state at all times. Thus, 
TBR may produce doubly-excited, and hence autoioniz- 
ing, Rydberg atoms with a considerably large rate. (This 
process is in close analogy to the production of autoion- 
izing Rydberg states by "isolated-core excitation" [HI-) 
For low enough principal quantum numbers, the autoion- 
ization rate of these states becomes comparable to and 
even exceeds the radiative decay rate of the excited core. 
For the case of strontium, the 1-averaged autoionization 
rate becomes important at n k 50. For electron tem- 
peratures of the present type of experiments, Rydberg 
atoms typically recombine into states with much higher 
n so that Auger ionization does not play a role initially. 
However, in the course of the gas evolution the Rydberg 
electron moves down the energy ladder by subsequent in- 
elastic electron-atom collisions. Since the energy shift of 
the core transition used for laser cooling is of the same 
order of magnitude as the autoionization rate |25| . the 
Rydberg atoms formed by TBR are still resonant with 
the cooling laser. Hence, even if the timescale of colli- 
sional deexcitation is longer than the lifetime of the core- 
excited state, the cooling laser will continue to drive the 
core transition so that the atom will be in the core-excited 
state for a significant fraction of time. Because the en- 
ergy connected with the core transition is of the order of 
10 4 K, each free electron produced by autoionization will 
rapidly leave the plasma volume. Hence, the combined 
action of the cooling laser, TBR and collisional deexci- 
tation is expected to remove electrons from the plasma 
and destroy the plasma, until recombination stops when 
the electron density has become too small. In order to 
suppress this electron loss one has to choose initial con- 
ditions which lead to a tolerable TBR rate ■ 

If time and length are measured in units of the initial 
ion plasma frequency lo p o and Wigner-Seitz radius do, 
respectively, the initial plasma state is characterized by 
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FIG. 2: Distribution of scaled inter-ionic distances after a 
time of t — 52 /is (io p ot = 240) compared to the calculated 
pair-correlation function of an OCP at T 4 = 700 |24|. The 
inset shows the ionic CCP calculated by different methods 
(see text for details). 



five parameters: the number of ions Ni, the initial elec- 
tron and ion CCPs, the value of T c corresponding to the 
Doppler temperature T c , and the ratio of the cooling rate 
to initial ionic plasma frequency lj p q. The exact values 
of T c and Ti(0) are not important for the plasma ex- 
pansion dynamics, since both are negligible compared to 
the electron temperature. According to Eq. for fixed 
P the time evolution of the scaled density only depends 
on the product of T e and N^ 3 (if Ti <C T e is neglected 
in the second equation of Eq. @), the expansion being 
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slower if Te-A^ becomes larger. Hence, in order to slow 
down the plasma expansion sufficiently for spatial cor- 
relations to develop, one may prefer to increase the ion 
number since reduction of the electron temperature be- 
comes ultimately incompatible with the objective of lim- 
iting TBR. A further constraint on the initial conditions 
arises from the fact that Ti « 10 3 in order to observe 
ordering effects. Typically, ion temperatures of the order 
of 1 mK can be achieved through Doppler cooling, mean- 
ing that the initial ion density must be at least about 10 s 
cm~ 3 . Finally, the cooling rate j3 must be of the order 
of 10 _1 cijpo to sufficiently slow down the plasma expan- 
sion that correlations can develop. Since [3 (oc 1/rrii) 
decreases faster than aj p0 (oc l/^/mi) with increasing ion 
mass, it is advantageous to consider relatively light ions, 
for which sufficiently high cooling rates can be obtained 
experimentally. We therefore choose Be ions for our sim- 
ulations, for which laser-cooling has been experimentally 
demonstrated earlier in nonneutral plasmas |2flj . 

The most striking result of our simulations is a lattice- 
type crystallization of the ions or even their arrange- 
ment in concentric shells if the plasma expansion is slow 
enough. The emergence of such an order depends sensi- 
tively on the initial conditions. We have found a lattice- 
type crystallization in the expansion of a plasma of 20000 
Be-ions with an initial density of 1.1 ■ 10 8 cm -3 and an 



electron temperature of T e = 29 K, cooled with a damp- 
ing rate of [3 = 0.15o-> p o to an ionic temperature of T c = 2 
mK. The value of Ti = 750 after 52/is, simply calculated 
from the ion temperature and the average Wigner-Seitz 
radius, suggests strong ordering of the ionic component. 
However, as pointed out in 9] the CCP calculated in 
this way may have no meaning as a measure of correla- 
tions, since the expanding plasma does not reach a global 
equilibrium. A more reliable quantity, namely the distri- 
bution of inter-ionic distances, is shown in Fig. [3 In 
order to account for the nonhomogeneity of the plasma 
we have scaled the inter-particle spacing by the Wigner- 
Seitz radius a/ oc determined by the local density between 
the corresponding particles. For comparison, the pair- 
correlation function obtained from the HNC-equations 
of a homogeneous one-component plasma is also 
shown. (However, the distribution function shown in 
Fig. should not be understood as a pair-correlation 
function in a strict sense, since in the present case the 
plasma is neither isotropic nor homogeneous.) The re- 
markable agreement shows that under present conditions 
the calculated CCP indeed indicates the degree of order 
in the expanding plasma and that the system has reached 
a state of local equilibrium far beyond the known crys- 
tallization limit of Ti w 174 0. In order to study the dy- 
namics of the crystallization process we have determined 
Ti also from the numerically obtained average correlation 
energy to geth er with an analytical approximation for this 
quantity [17|. A comparison of the CCPs calculated by 
both methods is shown in the inset of Fig. Initially 
there are large deviations between both calculations re- 
flecting the nonequilibrium character of the early plasma 
state. However, after some inverse plasma frequencies, 
the ion system reaches a local equilibrium and the CCPs 
obtained from the different methods become identical. 
At longer times both curves may diverge again due to a 
freezing out of ordered structures when the density be- 
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FIG. 3: Radial density after a time of t = 24 /is (uj p ot = 
110). The inset shows a two-dimensional cut through the 
plasma cloud, clearly revealing the formation of concentric 
shells. (For better contrast, cuts with x = 0, y — and z — 0, 
respectively, have been overlayed.) 
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FIG. 4: (Color online) The 5 th shell of Fig. 01 demonstrating 
significant intra-shell ordering. 

comes too small. 

If we further slow down the plasma expansion, by Ln- 

2/3 

creasing the product T e N i , the system exhibits for- 
mation of concentric shells rather than relaxation into 
lattice-type structures. This is demonstrated in Fig. [3J 
where the radial ion density is shown at t = 2A[is for a 
plasma with r e = 0.15 and N t — 50000 while all other 
parameters equal those used for Fig. [21 Besides this 
radial ordering into concentric shells, there are strong 
intra-shell correlations, also found in trapped nonneu- 
tral plasmas 0. In Fig. H we show one of the shells 
formed in the simulation of Fig.|3 The development of a 
hexagonal-lattice like ordering is evident, which is how- 
ever considerably disturbed by the curvature of the shell. 
A closer look on the emergence of the order shows that 
in the early stages of the plasma evolution a cubic-lattice 
like structure is formed. However, after some ten inverse 
plasma frequencies the ions rearrange to form concentric 
shells starting from the plasma center, where the den- 
sity is highest, in contrast to trapped nonneutral plas- 
mas where the shell formation was observed to proceed 
from the periphery to the center of the cloud . If the 
expansion is faster, as in the first example (Fig. [21, ion- 
ion collisions are less frequent and after the initial phase 
of local ordering the density is already too low for the 
rearrangement into shells, such that the lattice structure 
survives during the expansion. 

In summary, we have followed the long-time dynamics 
of laser-cooled, expanding ultracold plasmas on the ba- 
sis of a hybrid-MD simulation, allowing for a full treat- 
ment of the strongly coupled ion dynamics. The results 
show that cooling during the plasma expansion dras- 
tically modifies the expansion dynamics leading to an 
exotic type of plasma where the electron component is 
weakly coupled while the ion component shows strong 
coupling effects which manifest themselves in the devel- 
opment of lattice-like structures (short-range order) or 
even the formation of concentric shells (long-range or- 
der) depending on the expansion dynamics. Interesting 



questions concerning this novel system, like the behav- 
ior of ion collective modes or the influence of different 
density profiles on the details of the structure formation 
process, have to be addressed in future studies. 
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